home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / velovect.pro < prev    next >
Text File  |  1997-07-08  |  6KB  |  174 lines
  1. ; $Id: velovect.pro,v 1.12 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1983-1997, Research Systems, Inc.  All rights reserved.
  4. ;    Unauthorized reproduction prohibited.
  5.  
  6. PRO VELOVECT,U,V,X,Y, Missing = Missing, Length = length, Dots = dots,  $
  7.         Color=color, _EXTRA = extra
  8. ;
  9. ;+ 
  10. ; NAME:
  11. ;    VELOVECT
  12. ;
  13. ; PURPOSE:
  14. ;    Produce a two-dimensional velocity field plot.
  15. ;
  16. ;    A directed arrow is drawn at each point showing the direction and 
  17. ;    magnitude of the field.
  18. ;               
  19. ; CATEGORY:
  20. ;    Plotting, two-dimensional.
  21. ;
  22. ; CALLING SEQUENCE:
  23. ;    VELOVECT, U, V [, X, Y]
  24. ;
  25. ; INPUTS:
  26. ;    U:    The X component of the two-dimensional field.  
  27. ;        U must be a two-dimensional array.
  28. ;
  29. ;    V:    The Y component of the two dimensional field.  Y must have
  30. ;        the same dimensions as X.  The vector at point [i,j] has a 
  31. ;        magnitude of:
  32. ;
  33. ;            (U[i,j]^2 + V[i,j]^2)^0.5
  34. ;
  35. ;        and a direction of:
  36. ;
  37. ;            ATAN2(V[i,j],U[i,j]).
  38. ;
  39. ; OPTIONAL INPUT PARAMETERS:
  40. ;     X:    Optional abcissae values.  X must be a vector with a length 
  41. ;        equal to the first dimension of U and V.
  42. ;
  43. ;    Y:    Optional ordinate values.  Y must be a vector with a length
  44. ;        equal to the first dimension of U and V.
  45. ;
  46. ; KEYWORD INPUT PARAMETERS:
  47. ;      MISSING:    Missing data value.  Vectors with a LENGTH greater
  48. ;        than MISSING are ignored.
  49. ;
  50. ;    LENGTH:    Length factor.  The default of 1.0 makes the longest (U,V)
  51. ;        vector the length of a cell.
  52. ;
  53. ;    DOTS:    Set this keyword to 1 to place a dot at each missing point. 
  54. ;        Set this keyword to 0 or omit it to draw nothing for missing
  55. ;        points.  Has effect only if MISSING is specified.
  56. ;
  57. ;    COLOR:    The color index used for the plot.
  58. ;
  59. ;    Note:   All other keywords are passed directly to the PLOT procedure
  60. ;        and may be used to set option such as TITLE, POSITION, 
  61. ;        NOERASE, etc.
  62. ; OUTPUTS:
  63. ;    None.
  64. ;
  65. ; COMMON BLOCKS:
  66. ;    None.
  67. ;
  68. ; SIDE EFFECTS:
  69. ;    Plotting on the selected device is performed.  System
  70. ;    variables concerning plotting are changed.
  71. ;
  72. ; RESTRICTIONS:
  73. ;    None.
  74. ;
  75. ; PROCEDURE:
  76. ;    Straightforward.  Unrecognized keywords are passed to the PLOT
  77. ;    procedure.  
  78. ;
  79. ; MODIFICATION HISTORY:
  80. ;    DMS, RSI, Oct., 1983.
  81. ;    For Sun, DMS, RSI, April, 1989.
  82. ;    Added TITLE, Oct, 1990.
  83. ;    Added POSITION, NOERASE, COLOR, Feb 91, RES.
  84. ;    August, 1993.  Vince Patrick, Adv. Visualization Lab, U. of Maryland,
  85. ;        fixed errors in math.
  86. ;    August, 1993. DMS, Added _EXTRA keyword inheritance.
  87. ;    January, 1994, KDB. Fixed integer math which produced 0 and caused
  88. ;                    divide by zero errors.
  89. ;    December, 1994, MWR. Added _EXTRA inheritance for PLOTS and OPLOT.
  90. ;    June, 1995, MWR. Removed _EXTRA inheritance for PLOTS and changed
  91. ;             OPLOT to PLOTS.
  92. ;       September, 1996, GGS. Changed denominator of x_step and y_step vars. 
  93. ;-
  94. ;
  95.         on_error,2                      ;Return to caller if an error occurs
  96.         s = size(u)
  97.         t = size(v)
  98.         if s[0] ne 2 then begin 
  99. baduv:   message, 'U and V parameters must be 2D and same size.'
  100.                 endif
  101.         if total(abs(s[0:2]-t[0:2])) ne 0 then goto,baduv
  102. ;
  103.         if n_params(0) lt 3 then x = findgen(s[1]) else $
  104.                 if n_elements(x) ne s[1] then begin
  105. badxy:                  message, 'X and Y arrays have incorrect size.'
  106.                         endif
  107.         if n_params(1) lt 4 then y = findgen(s[2]) else $
  108.                 if n_elements(y) ne s[2] then goto,badxy
  109. ;
  110.         if n_elements(missing) le 0 then missing = 1.0e30
  111.         if n_elements(length) le 0 then length = 1.0
  112.  
  113.         mag = sqrt(u^2.+v^2.)             ;magnitude.
  114.                 ;Subscripts of good elements
  115.         nbad = 0                        ;# of missing points
  116.         if n_elements(missing) gt 0 then begin
  117.                 good = where(mag lt missing) 
  118.                 if keyword_set(dots) then bad = where(mag ge missing, nbad)
  119.         endif else begin
  120.                 good = lindgen(n_elements(mag))
  121.         endelse
  122.  
  123.         ugood = u[good]
  124.         vgood = v[good]
  125.         x0 = min(x)                     ;get scaling
  126.         x1 = max(x)
  127.         y0 = min(y)
  128.         y1 = max(y)
  129.     x_step=(x1-x0)/(s[1]-1.0)   ; Convert to float. Integer math
  130.     y_step=(y1-y0)/(s[2]-1.0)   ; could result in divide by 0
  131.  
  132.     maxmag=max([max(abs(ugood/x_step)),max(abs(vgood/y_step))])
  133.     sina = length * (ugood/maxmag)
  134.     cosa = length * (vgood/maxmag)
  135. ;
  136.         if n_elements(title) le 0 then title = ''
  137.         ;--------------  plot to get axes  ---------------
  138.         if n_elements(color) eq 0 then color = !p.color
  139.         x_b0=x0-x_step
  140.     x_b1=x1+x_step
  141.     y_b0=y0-y_step
  142.     y_b1=y1+y_step
  143.         if n_elements(position) eq 0 then begin
  144.           plot,[x_b0,x_b1],[y_b1,y_b0],/nodata,/xst,/yst, $
  145.             color=color, _EXTRA = extra
  146.         endif else begin
  147.           plot,[x_b0,x_b1],[y_b1,y_b0],/nodata,/xst,/yst, $
  148.             color=color, _EXTRA = extra
  149.         endelse
  150. ;
  151.         r = .3                          ;len of arrow head
  152.         angle = 22.5 * !dtor            ;Angle of arrowhead
  153.         st = r * sin(angle)             ;sin 22.5 degs * length of head
  154.         ct = r * cos(angle)
  155. ;
  156.         for i=0,n_elements(good)-1 do begin     ;Each point
  157.                 x0 = x[good[i] mod s[1]]        ;get coords of start & end
  158.                 dx = sina[i]
  159.                 x1 = x0 + dx
  160.                 y0 = y[good[i] / s[1]]
  161.                 dy = cosa[i]
  162.                 y1 = y0 + dy
  163.         xd=x_step
  164.         yd=y_step
  165.                 plots,[x0,x1,x1-(ct*dx/xd-st*dy/yd)*xd, $
  166.             x1,x1-(ct*dx/xd+st*dy/yd)*xd], $
  167.                       [y0,y1,y1-(ct*dy/yd+st*dx/xd)*yd, $
  168.             y1,y1-(ct*dy/yd-st*dx/xd)*yd], $
  169.                       color=color
  170.                 endfor
  171.         if nbad gt 0 then $             ;Dots for missing?
  172.                 PLOTS, x[bad mod s[1]], y[bad / s[1]], psym=3, color=color
  173. end
  174.